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ABSTRACT 



f 



Aims. The aim of this paper is to investigate the properties of accretion disks threaded by a weak vertical magnetic 
field, with a particular focus on the interplay between MHD turbulence driven by the magnetorotational instability 
(MRI) and outflows that might be launched from the disk. 

Methods. For that purpose, we use a set of numerical simulations performed with the MHD code RAMSES in the 
framework of the shearing box model. We concentrate on the case of a rather weak vertical magnetic field such that 
the initial ratio /3o of the thermal and magnetic pressures in the disk midplane equals 10*. 

Results. As reported recently, we find that MHD turbulence drives an efficient outflow out of the computational box. We 
demonstrate a strong sensitivity of that result to the box size: enlargements in the radial and vertical directions lead to 
a reduction of up to an order of magnitude in the mass-loss rate. Such a dependence prevents any realistic estimates of 
disk mass-loss rates being derived using shearing-box simulations. We find however that the fiow morphology is robust 
and independent of the numerical details of the simulations. Its properties display some features and approximate 
invariants that are reminiscent of the Blandford & Payne launching mechanism, but differences exist. For the magnetic 
field strength considered in this paper, we also find that angular momentum transport is most likely dominated by MHD 
turbulence, the saturation of which scales with the magnetic Prandtl number, the ratio of viscosity and resistivity, in 
a way that is in good agreement with expectations based on unstratified simulations. 

Conclusions. This paper thus demonstrates for the first time that accretion disks can simultaneously exhibit MRI-driven 
MHD turbulence along with magneto-centrifugally accelerated outfiows. However, in contradiction with previously 
published results, such outflows probably have little impact on the disk dynamics. 



1. Introduction 

In the last two decades, MHD turbul ence mediated by 
the magnet orotational instability (MRI; iBalbus fc HawlevI 
Il99ll Il998f) has been extensively studied as the most 
likely mechanism responsible for angular momentum 
transport in accretion disks. Yet the natural and simple 
configuration of an initially uniform vertical magnetic 
field threading an isothermal vertically stratified disk has 
been relatively neglected compared to the large number of 
papers that have been published reporting studies of other 
configurations, such as unstratified disks, or stratified disks 
with a toroidal magnetic field or without any net flux 



their numerical codes or the catastrophic disruption of the 
disk in a few dyna mical times. T h ese di sturbing findings, 
later confirmed by I Turner et al.l (|2006[ ). led researchers 
to focus on other, more user-friendly, configurations. 
This is unfortunate because a configuration with a net 
poloidal magnetic field is the most natural to explain 
the widely observed o ccurrence of iets and win ds in 
astrophysical systems (jSpruitl Il996l : IWardld Il997l and 
references therein). Magnetic fields are also believed to 
be of 



(IStone et al.l l996': Brand enburg 



Zieglcr fc Ritdigeij |2000F 



Hirose et al.1 l2006t 



Oishi &: Mac Low! 12009): IShi et all l2010t iDavis et al 
iSi 



Johansen &: Levir] 
«^ 

I; IF 



cess (iHennebelle & Fromane 20081: Hennebelle & Ciardil 


20091: Commercon 


et al.l 120101: IHennebelle et al.l 120111: 


Seifried et al.l 2011 


; iJoos et al.l|2012f). Even if diffusion is 



19£5[^ Miller fc Stone required to remove much of the magnetic fiux during that 



Fleming fc Stone 
20n8t 



process (iMellon fc Lill2 009l: 'Krasnopols ky et al.l[20T 
iLi et al.ll2011l: ISantos-L ima et al. 2012: iDapp et al 



I 



ionl; 

mi, 



20101: iFlaig eralTl2010l: Tsimon et al. l2011bl ial:lFlaig et al. 
20121) . This is largely the result of computational diffi- 
cul ties. Indeed, early simula tions performed in th e 1990s 
by IStone et al.l ([l996h and iMiller fc Sto^ (I2OOOI) found 
that the configuration of a net vertical magnetic field 
in a stratified disk resulted in either the breakdown of 



it is plausible that some net poloidal flux is retained and 
carried down to the scale of the disk during its formation. 
The resulting strength of that potential magnetic fleld is 
highly uncertain, but might be in the approximate range 
10"^ to a few Gauss (Wardlc 2007). This translates into /3 
values (the ratio between thermal and magnetic pressure) 
ranging from 10^ to 1. 
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Recently, however, the net vertical field configuration 
has received renewed interest. Using local simulations per- 
form ed in the shearing-box ap proximation in that situa- 
tion, ISuzuki fc Inutsukal (|200a hereafter SI09) found that 
accretion disks can drive winds able to remove significant 
amounts of mass within dynamical timescales. By varying 
the initial value of the plasma /? parameter, the ratio of mid- 
plane thermal and magnetic pressures, between 10*^ and in- 
finity (i.e. no vertical magnetic field), they found that the 
mass outflow rate driven by such winds is an increasing 
function of the magnetic field strength. For the strongest 
field they investigated, /3 = 10''^, the outflow is so powerful 
that the disk is emptied in less than 1 00 local orbits! As 
investigated later by the same authors (jSuzuki et al.ll2010l 
hereafter SMIlO), such strong disk outflows have important 
consequences for disk evaporation and could also help halt- 
ing inward planet migration. This is an important result 
that, as such, requires confirmation and follow-up. This is 
the first goal of this paper. 

Another important unresolved issue under these condi- 
tions is the saturation of the turbulence its elf. Using un- 
stratified shearing boxes, early simulations (jHawlev et al.l 
soon established that the rate of angular momen- 
tum transport scales like the initial vertical magnetic 
field strength. Note however that this s caling has recently 
been questione d by various authors (jBodo et al.l 1201 it 
IUzdensk"v 20121). In addition, extensive parameter su rveys 
( Lesur fc Longarettil 120071 : iLongaretti fc Lesurl |2010( ) per- 
formed using the same setup also demonstrated a strong 
dependence of the transport with the magnetic Prandtl 
number Pm, the ratio of kinematic viscosity to ohmic resis- 
tivity. How these results arc modified when taking vertical 
stratification into account is still unknown. Performing a 
first step in that direction is the sec o nd go al of that pa- 
per. As recently argued by lUzdenskvl (j2012|) . the question 
of the saturation of MHD turbulence in such a situation is 
not only important for angular momentum transport itself, 
but is also linked to issues such as disk coronal heating and 
magnetic fiux transport. 

The plan of the paper is as follows. In section 12.11 we 
present the numerical method we used and describe the 
extensions of the basic scheme that were required for the 
simulations to run without failure. The immediate result of 
these simulations is a confirmation of the results reported 
by SI09: an MHD turbulent disk threaded by a net vertical 
flux drives outflows that efficiently remove mass from the 
computational box. In section [3] we focus on the sensitivity 
of the results to details of the numerical setup, such as 
resolution and box size. Finally, in section 21 we make a 
connection with existing disk wind theories and study the 
saturation properties of the turbulence, before discussing 
some aspects of our results (section [5]) and perspectives for 
future work (section [5]). 



2. Methods 

2.1. Numerical setup 

The numerical simulations we present in this paper 
use a setup similar to that of IStone et"al] ()1996t ) and 
iBrandenburg et all (|l995f ). The MHD equations are solved 
in a Cartesian coordinate system (x, y, z) with unit vec- 
tors {i,j,k) rotating with angular velocity il around a 
central mass. This is the classical shearing-box model 



(iGoldreich fc Lvnden-Belll 119651: iHawlev et all I1995D . As 

described by previous authors, these equations include a 
vertical component of the gravitational force and possibly 
also viscous and ohmic dissipation (see also iDavis et al.l 
I2OIOD . For simplicity, we consider only the case of an 
isothermal gas with sound speed cq. 

As described in the introduction, the aim of the present 
paper is to examine configurations in which the vertical 
magnetic flux (conserved during the simulation because of 
the shearing box symmetries) does not vanish. This is done 
by initializing the magnetic field as a uniform vertical field 
JB ~ Bok whose strength is parametrized by the midplane 
plasma parameter /3 defined according to 



(1) 



where po is the value of the gas density p in the disk 
midplane, and we work in electromagnetic units such that 
Po = 1. In this paper, we consider the case /3o = 10^. This 
purely vertical magnetic field is superposed on the initial 
hydrostatic disk configuration at the beginning of the sim- 
ulations, along with velocity fluctuations that trigger the 
growth of the MRI. 

We used standar d shearing-box boun dary conditions in 
the radial direction (jHawlev et al.l 119951 ) and periodic con- 
ditions in the azimuthal direction. At the vertical bound- 
aries, however, we used modified outflow boundary condi- 
tions. The density is extrapolated assuming vertical hydro- 
static equilibrium. Zero-gradient boundary conditions are 
applied on horizontal velocities and on the vertical momen- 
tum, when matter is outflowing (otherwise, the vertical ve- 
locity is set to zero in the ghost cells). Finally, the magnetic 
field is forced to be vertical in the ghost cells. 

The set of equations just described are solved using a 
version of the code RAMSES (|Tevssiedl2002l:lFromang et al.l 
I2OO6D that solves the MHD equations on a uniform grid 
in the shea r ing b ox. It was quantitatively tested by 
iLatter et al.l (l2010l ) and the results obtained in 3D sim- 
ulations of MRI-driven turbulence were shown to com- 
pare successfully w ith those obtained with the code Athena 
(|Stone et al.l I2008D by iFromang fc Stond (|2009t ). In the 
present paper, several simulations are performed in a ra- 
dially extende d box (larger than the disk scale height H). 
iJohnson et al. (2008) have shown that position-dependent 
dissipation introduces artificial radial variations of the den- 
sity and the stress in this case. As shown by the same au- 
thors, this effect can be eliminated by using an M HD ex- 
tensi on of the FARGO orbital advection algorithm (jMassetl 
I2000h . For the purpose of this work, we have implement 
su ch an extension, fol l owing the method recently described 
bv lStone fc Gardineil (|2010t ) for finite- volume methods. 

2.2. Extensions of tlie basic scheme 

Shearing-box numerical simulations that start with vertical 
magnetic flux lead to vigorous turbulence. Strong magnetic 
flelds develops in the upper layers of the disk and can cause 
the code to crash randomly during a run. For this reason, 
a number of modifications had to be implemented for the 
simulations to proceed robustly. We describe them here. 

First, RAMSES uses the MUSCL- Hancock algorithm to 
integrate the MHD equations (jTorol 119971 ). Spatial slopes 
of all variables are required for the scheme to be of sec- 
ond order in space and time. Simulations performed with 
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RAM SES normal ly use the MinMod or MonCen slope lim- 
iters (iTorol 1199711 . However, we found that both limiters 
systematically produce evacuated regions where the Alfvcn 
velocity is enormous, resulting in extremely small time steps 
that cause the code to halt. We therefore use d the multi- 
dimensional limiter described by ISureshI ()2000t ) that is able 
to avoid such problems. 

However, this turned out to be insufficient in radially 
extended boxes. Two additional modifications were imple- 
mented. First, at those cell interfaces (usually located in 
the disk corona) where the magnetic pressure exceeds the 
thermal pressure by more than three orders of magnitude, 
we used the Lax -Friedrichs Riemann solv er. The HLLD 
Riemann solver (iMivoshi fc Kusanol l2005l) is used at all 
other locations. In addition, a mass diffusion source term 
was added to the continuity equation. While conserving 
mass, it helps fill the density holes described above and 
prevent large Alfven speeds from appearing. We used the 
same analytical exp ression for the diffusion coefficient as 



same analytical exp i 
iGressel et al] (poTl : 



D = Do 



Po 



1/4 



(2) 



with Do = 10~^c§/ri and Cdj,„ ~ 5. Such parameters yield 
a grid Reynolds number for mass diffusion of order unity in 
the disk upper layers. 

When used simultaneously, the modifications to the 
standard scheme described in this section permit a robust 
execution of the code until completion of the simulations. 
The modifications they induce of the properties of the flow 
are undetectable. 



2.3. Run properties 

Wc summarize in Table [1] the properties of the simulations 
presented in this paper. Table [1] shows the model label 
(column 1), the box size (column 2), the grid resolution 
(column 3) and the run duration (column 4). Columns 5 
and 6 report, where applicable, the Reynolds and magnetic 
Reynolds numbers we used. They are respectively defined 
by: 

,?e=£2£andi?™=^, 
v rj 

where H = co/fl is the disk scale height while v and rj 
are the kinematic viscosity and ohmic resistivity (magnetic 
diffusivity) respectively. Finally, the average dimensionless 
mass-loss rate per unit area through the upper and lower 
surfaces of the box is reported in column 7. It is calcu- 
lated according to 



1 {pVz)top 

2 



{Pvz) 



hot 



PoCoLxLy 



(3) 



where {pVz)top and {pVz)bot are respectively the horizontally 
integrated mass fluxes through the lower and upper surfaces 
of the box. 



3. Numerical issues 

In all the simulations we performed using the setup de- 
scribed above, we found that strong disk outflows similar to 
those described by SI09 and SMIlO develop. This indicates 



that outflows are a robust outcomes of MHD shearing-box 
simulations of stratifled disks threaded by a vertical mag- 
netic fleld. In this section we analyse in detail the influ- 
ence of various computational parameters, such as vertical 
boundary conditions, resolution and box size, on the prop- 
erties of these outflows. 



3.1. Comparison with \Suzuki Si Inutsukd (200^) 

As described in section [2.11 we use vertical boundary con- 
ditions that are different from SI09. Indeed, these authors 
prescribed outflow conditions by adopting only outgoing 
MHD ch aracteristics at their top and bottom boundaries 
(jSuzuki fc Inutsuka 2006). Our choice of different boundary 
conditions provides an opportunity to test the sensitivity of 
the outflows to this aspect of the numerical setup. To assess 
this possibility in a quantitative way, we reproduce their 
model as closely as possible by performing a simulation in 
a box of the same size {Lx,Ly, L^) = {V^H, Ay/^H, 8V2H) 
and with identical resolution {Nx,Ny,Nz) = (32,64,256). 
The factor in the box size appears because of our differ- 
ent definition of H compared to SI09 (we will later return to 
box sizes that are integer numbers of H) . This model is la- 
belled SIOQrun (see Table [T]) . Note that the simulation was 
only evolved for about 20 orbits, which is short by modern 
standards. However, 20% of the disk mass has escaped the 
computational box by that time, and we found that this is 
large enough to affect the subsequent flow properties. The 
mass-loss rate per unit area rhw we measured in this model, 
averaged in time between 5 and 20 orbits, is 2.4 x 10"'^. SI09 
do not quote explicitly their values, which can however be 
read off from their figures. In their figure 5, SI09 report 
2mu, ^ 5-6 X lO^'^. Thus there is good agreement with 
our results, which suggests the particular choice of bound- 
ary conditions has only a limited effect on the quantitative 
properties of the fiow. 



3.2. Influence of the vertical resolution 

A second peculiar a spect of SI09 is th e low r esolution they 
used, as noted by iGuan fc Gammid ()2011|) . Indeed, the 
most unstable MRI mode has a vertical wave number given 
by k^VAz ~ in an unstratified box. For (3 = 10^, as con- 
sidered in the present paper and by SI09, this corresponds 
to a wavelength of about one tenth of the disk scale height. 
This wavelength is not adequately resolved in model ST09a 
in which we use 32 cells per H. In stratified disks, the eigen- 
mode structure is more co mplex and cannot b e reduced to 
a single wave number, but iLatter et al.l (j2010t) have shown 
that about 200 cells per scale height are required in this case 
to properly resolve the fastest growing eigcnmode across the 
whole disk. In order to assess the importance of the verti- 
cal resolution on the existence and strength of the outffows 
that develop in model ST09a, we computed four models 
{Ideal256 to Ideal2048), gradually increasing the vertical 
resolution from = 256 to Nz = 2048 while keeping 
all other parameters fixed. In figure [TJ we plot the time- 
evolution of the total gas mass M in the computational 
box, normalized by its initial value Mq. We find it is simi- 
lar in the four models, regardless of the number of vertical 
cells. The values of rriw for the four models are reported in 
Table [TJ They confirm the results of figure [TJ as all values 
agree within ±10%. It would be tempting to associate this 
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Model 


Box size 


Resolution 


Run time (orbits) 


Re 


Rm 






ZRmSOOO 


{4:H,4H,H) 


(256, 128, 64) 


100 


3000 


3000 


_ 




ZRmlSOO 


[AH,AH,H) 


(256, 128, 64) 


100 


3000 


1500 


- 




SI09run 


{V2H, 4:V2H, 8V2H) 


(32,64,256) 


20 


_ 


_ 


2.4 X 10" 


-3 


Ideal256 


{H, AH, lOH) 


(32, 64, 256) 


20 






2.6 X 10" 


-3 


Ideal512 


{H, AH, IQH) 


(32,64,512) 


20 






2.7 X 10" 


-3 


Ideall024 


[h, ah, IQH) 


(32, 64, 1024) 


20 


_ 


_ 


2.3 X 10" 


-3 


Ideal2048 


{H, AH, WH) 


(32, 64, 2048) 


20 


_ 


_ 


2.8 X 10" 


-3 


DiffulH 


[H, AH, IQH) 


(64, 128, 640) 


30 


3000 


3000 


2.2 X 10" 


-3 


Diffu2H 


{2H, AH, IQH) 


(128, 128,640) 


30 


3000 


3000 


1.5 X 10" 


-3 


Diffu4H 


{AH, AH, IQH) 


(256,128,640) 


40 


3000 


3000 


1.1 X 10" 


-3 


DiffuSH 


{8H, 8H, IQH) 


(512,256,640) 


10 


3000 


3000 


1.1 X 10" 


-3 


Tall4H 


{AH, AH, 2QH) 


(256, 128, 1280) 


20 


3000 


3000 


3.2 X 10" 


-4 


Diffu4Ha 


{AH, AH, IQH) 


(256, 128,640) 


40 


3000 


1500 


1.2 X 10" 


-3 



Table 1. Properties of the runs described in this paper. The first column gives the model label, while the next two 
columns give the size of the computational domain {L^, Ly, L^) and the resolution {N^, Ny, Nz) of that run. We give 
in column four the duration (in orbits) of each simulation. Dissipation coefheients, expressed in terms of the Reynolds 
and magnetic Reynolds numbers, are given in the following two columns. Finally, the last column gives, in dimensionless 
form, the rate m^j at which matters escapes the computational domain per unit surface area. 



1.05 
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10 

Time (in orbits) 



Fig. 1. Time evolution of the total relative mass in the box 
when the number of cells in the vertical direction is varied; 
Nz = 256 {blue line), = 512 {green line), = 1024 
{red line) and = 2048 {light blue line). All other run 
parameters arc fixed. 



result with the properties of the eigenmodes along with the 
fact that disk outflows are launched at z '--^ zL2H, as done 
by SI09 and SMIIO. Indeed, MRI modes have a shorter 
wavelength nea r the midplane an d a longer wavelength at 
greater height (jLatter et al.ll2010[) . Since outflow rates are 
not set by conditions at the midplane, the mode structure 
at the mid plane may be under resolved with little impact. 
In addition, iLatter et al.l ()2010[) have also shown that those 
modes are not only solution of the linear equation, but also 
nonlinear solutions in the limit of large /3, such that some of 
their properties might affect the nonlinear flow dynamics. 
This reasoning would also agree with the claim of SI09 that 
such disk outflows are due to the breakup of channel modes 
in the disk corona. It is indeed likely that the properties of 
those channel with help the build up of strong field in the 



disk atmosphere during the linear stage of the MRI. That 
having being said, the dynamics of the disk atmosphere 
during the remainder of the simulations is highly nonlinear 
and fluctuating. It might share some properties of the dy- 
namics of channel modes, but a direct relationship between 
such outflows and linear MRI modes remains difficult to 
assess with certainty. 

3.3. Explicit dissipation coefficients 

Even though the previous section has established that the 
mass outflow rate is independent of the vertical resolu- 
tion, these simulations, in order to stabilize the numer- 
ical scheme, still rely on grid dissipation, the form and 
effects of which are more difficult to quantify than those 
of explicit dissipation. As explained in the introduction, 
it is also well known that the saturation properties of 
the tu rbulence depend strongly on the dissipatio n coeffi- 
cients (|Lesur fc Longa retti 2007: Longaretti fc Le sur 203; 
'Fromang et al."2007'; Fromang 20l0':'Si mon &: Hawlev)l200a 
Davis et aL 2010; Simon et al. 2011b. ), regardless of the 
magnetic field configuration. 

The introduction of explicit dissipation coefficients also 
makes the problem mathematically well posed. In contrast 
to the low-vertical-resolution simulations described above 
(for which numerical dissipation influences the linear phase 
of the instability), when viscosity and Ohmic resistivity 
are included the calculated growth rates and morphology 
of MRI normal modes can be accurately reproduced in 
the simulations, providing a valuable code test. In the ap- 
pendix, we describe the calculation of these modes as a 
function of Re an d Rm, extending the method used by 
ILatter et all (|2010f ) in the ideal MHD limit. As shown in 
fig- lA.li the linear modes display large-scale oscillations at 
z ~ zt3H regardless of Re and Rm, while the small-scale 
oscillations near the midplane are very much influenced by 
the magnitude of the dissipation coefficients. This structure 
could explain the relative independence of the flow struc- 
ture on resolution described in section [3.21 

One obvious limitation associated with explicit dissipa- 
tion coefficients, though, is the tremendous computational 
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1.05 




Time (in orbits) 

Fig. 2. Time-evolution of the total relative mass in the box 
when the radial extent of the box is varied from ^ H 
{blue line) to = 2H (green line) and = 4ff (red line). 



costs associated with the grid resolution that is required 
for the dissipative lengthscales to be properly resolved. To 
keep the resolution requirements reasonable, we adopt in 
the remainder of this paper rather modest values for the 
dissipation coefficients such that Re = 3000 for all simula- 
tions. We consider two values for Rm, namely Rm = 3000 
in our fiducial model and Rm = 1500 in section [42] Past 
experi ence gained from uns tratified simulations (see for ex- 
ample |^o^Sn^|er^L|[2003) suggests that 64 cells per disk 
scale height are sufficient to capture correctly the flow prop- 
erties at all relevant scales. Before moving to a systematic 
study of the flow in such a viscous and resistive model, we 
briefly compare the outflow rate obtained in one such case 
(which wc label model DiffulH for future reference) with 
the cases without explicit dissipation described above, all 
other parameters being unchanged. As shown in Table [1] 
TTiu, is similar in all cases, regardless of the nature of the 
dissipation. This is in agreement with the above result that 
mass outflow rates are independent of the resolution, when 
explicit dissipation is absent. We will therefore adopt the 
setup of model DiffulH as our fiducial model in the fol- 
lowing and we now move to a extensive study of the disk 
outflow properties and saturation level of the turbulence 
using these parameters. 

3.4. Influence of box size 
3.4.1. Radial box size 

In a suite of numerical simulations of unstratifled shea ring 
boxes threaded by a net vertical flux. lBodo et al.l(|2008^ dis- 
covered that the radial extension of the domain can signifl- 
cantly affect the properties of the turbulence. They demon- 
strated that channel flows develop in boxes of unit radial 
extent {L^ = H) while their influence is reduced in boxes 
having larger L^. Since SI09 argue that the breakup of these 
channel flows drives structured disk winds, it is natural to 
ask whether the radial extent of the shearing box affects the 
disk outflow properties described above. This is the purpose 
of the present section. 





A 

T \ 






6 -4 


-2 2 4 f 
Z/H 



Fig. 3. Vertical proflle of the ratio /3 between the horizon- 
tally averaged thermal pressure and the horizontally aver- 
aged radial pressure. The blue curve corresponds to model 
DiffulH and results from a further time average between 
t = 5 and t = 15 orbits. The red curve was obtained aver- 
aging the data of model Diffu4Hhetween t = 10 and t ~ 30. 
For each model, the vertical arrows report the approximate 
wind launching location where (3 = 1 (SI09). 



To do so, we performed a series of simulations, succes- 
sively doubling the radial extent of the box from = IH 
(model DiffulH) to = 4H (model Diffu4H). For each 
model, the time-evolution of the total mass inside the com- 
putational domain is shown in figure O The mass- loss rate 
gradually decreases as is increased. This is confirmed 
by Table [U where rhw is found to be reduced by a fac- 
tor of 2 in the largest model. Such a lower mass-loss rate, 
apart from changing the properties of the outflow itself, 
also allows the simulation to be run for longer (and more 
statistics of the turbulence to be accumulated) before the 
reduction of the disk surface density starts affecting the 
disk structure. This is the reason why we could run model 
Diffu4H for 40 orbits. A model with = 8H (model 
DiffuSH) was also run and suggested numerical conver- 
gence for larger radial extent. However, its large resolution, 
{N^,Ny,N^) = (512,256,640), resulted in huge computa- 
tional requirements that prevented the simulation from be- 
ing run for more than 10 orbits. Clearly, this question needs 
further investigation but our preli minary result i s not in dis- 
agreement with the conclusions of iBodo et al.l (j2008D . who 
find good agreement between their = 4: and = S 
runs. 

For a better understanding of the differences between 
these simulations, we turn to an examination of the wind 
launching location. SI09 have shown that this position can 
be evaluated to a good approximation by computing the lo- 
cation where the ratio /3 between the horizontally averaged 
thermal pressure and the horizontally averaged magnetic 
pressure falls below unity. The vertical profile of /3 is plot- 
ted in figure [3] for models DiffulH and Diffu4H. For both 
models, /? is large in the turbulent midplane and decreases 
upwards. The location where it falls below unity is re- 
ported using vertical arrows. For the smallest radial box 
size, we measured = ~2.1H and 1.4_ff respectively be- 
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low and above the equatorial plane. The same values are 
Zu] = —2.3H and 2.6H when = 4:H. Thus, in bigger 
boxes, disk outflows arc launched at higher altitude, where 
the gas density is smaller, resulting in a smaller mass-loss 
rates. Ultimately, the higher location of the launching re- 
gion can be be traced back to the toroidal magnetic field 
being weaker in the disk upper layers in bigger boxes, as 
shown by the larger 13 values obtained in figure [3] for model 
Diffu4H. The physical mechanism responsible for this re- 
duction is uncertain, but could be due to t he weaken- 
ing of coherent structures as demonstrated by iBodo et al] 
(|2008[ ) in unstratified boxes. This weakening indeed results 
in more disordered magnetic fields (seen in the simulation 
as a larger ratio between the fluctuating and means parts 
of the azimuthal magnetic field) that reduce the efficiency 
with which outflows are accelerated. 

3.4.2. Vertical box size 

In numerical simulations, it is desirable for the outflow ve- 
locity to exceed all wave velocities before leaving the com- 
putational domain. When this occurs, no signal can prop- 
agate downwind, which guarantees that the results are in- 
dependent of the box size and the boundary conditions^ 
In the case of the magnetized outflows studied here, the 
relevant wave velocities are those of the Alfven wave and 
the slow and fast magnetosonic w aves, each of which gives 
rise to a critical point in the flow ()Spruitlll996l ). Given the 
complexity of the flow, we will not consider wave velocities 
parallel to the fluid streamlines in this paper, but rather will 
focus on vertical velocity and wave velocities relative to the 
fluid. In the published simulations of SI09 and SMIlO, no 
information is given as to whether the critical points are 
crossed, but this is probably not the case, because the flow 
velocity remains subsonic in some of the cases they report 
(see for example figure 3 of SI09). 

In the present section, we investigate this issue by com- 
paring the results of model Diffu4-H, described above, with 
those of model TaU4H. The latter differs from the for- 
mer only in its vertical extent, namely Lz = 20H. We 
emphasize that the resolution of that simulation is large: 
{Na;,Ny,N:,) = (256,128,1280), which amounts to 40 mil- 
lion cells. Thus the computational requirements are signifi- 
cant and that model was evolved for only 20 orbits. This is 
however sufficient to derive a reliable estimate of the mass- 
loss rate, as shown by both figures [1] and [5] for the models 
we have described so far. We found that the mass outflow 
is significantly reduced when doubling the vertical extent 
of the box. Indeed, table [T] shows that ifiw is reduced by a 
factor of 3.4. Altogether, this means that while 20% of the 
mass in the box is lost within 20 orbital times in the sim- 
ulations of SI09, it would take about 140 orbits to lose the 
same amount of mass in model Tall^H. The consequences 
for disk evolution would surely be significant. 

This dependence on box size and its relation to the crit- 
ical points is further confirmed in figure 21 where the verti- 
cal profile of the vertical velocity (normalized by the sound 
speed) is compared with the vertical profile of the verti- 
cally propagating wave velocities (also normalized by cq) 



^ We note, however, that such a situation does not guarantee a 
successful launching of the outflow, as this depends on its ability 
to escape from the gravitational potential. We will come back 
to this point in section [5^ 




z/H 

Fig. 4. Vertical profiles of the vertical fluid and wave ve- 
locities. Horizontally and time-averaged fluid velocity 
(blue curve), slow magnetosonic speed v^z {light blue dashed 
curve), Alfven speed vaz {green dashed curve) and fast 
magnetosonic speed Vfz {red dashed curve), in units of cq. 
The top panel shows results obtained in model DiffuJ^H 
{Lx — iH, Lz ~ lOH) while the bottom panel is for the 
model TalUH {L.^ ^AH,Lz^ 2QH). 



for models Diffu^H {top panel) and Tall^H {bottom panel). 
The latt er are given b y the solutions of the following ex- 
pression (I0gilviell2012l ): 

[v'-vlz]y~{cl + vl)v'^clvlz]=0, (4) 

where Va ~ B/ ^ is the Alfven velocity. Two comments 
can be made. In both panels, it is seen that only the slow 
and Alfven points are crossed, while the fast magnetosonic 
point is, suspiciously, reached at the vertical boundary it- 
self. The fact that this is the case in both models suggest 
that this is no coincidence, but is rather a numerical ar- 
tifact probably associated with the shearing-box model or 
simply with the limited vertical extent of the computational 
domain. We note that such a failure to cross the fast mag- 
netosonic point is not unlike published simul ations of mag- 
netically driven jets (jCasse fc Ferreirall2000l see their fig. 
1). A second remark is that the slow and Alfven points lie 
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systematically higher up in the disk atmosphere in the taller 
model. The non-convergence of the location of those points 
might explain the dif ferent mass outflow rates measured in 
the two simulations (jSpruitllf 9961 ). 

4. Flow properties 



As shown in both sections 13.4.11 and 13.4.21 it appears that 
many details of the computational setup affect the proper- 
ties of outflows from MHD turbulent disks. The disappoint- 
ing but unavoidable conclusion from these results is that 
the disk mass-loss rate due to such outflows cannot be ro- 
bustly inferred in the framework of the shearing-box model. 
Nevertheless, we found that a number of properties of the 
flow are displayed by all simulations. The purpose of the 
present section is to describe these features and to analyse 
quantitatively the angular momentum transport associated 
with the turbulence. 

However, one should be cautious when doing so. As al- 
ready mentioned, the strong mass loss due to the outflow 
changes the disk structure itself on timescales of a few tens 
of orbits. For example, in model DiffulH, we discovered 
that a strong mean toroidal field, with a strength close to 
equipartition, emerged in the disk midplane after 30 orbits, 
perturbing its structure thereafter. No such perturbation 
was evident in model Diffu4H. In fact, all the diagnostics 
we considered suggest the disk is in a quasi-steady state 
during the first 20 to 25 orbits that follow the establishment 
of fully developed turbulence. This is why we will restrict 
our analysis in this section to the model having = 4i/. 
Unless otherwise specified, time averages of the data will 
cover the range 5 < i < 30. 

4.1. Disk atmosphere: outflow structure 

We first illustrate the flow structure with the help of flg- 
ure O Four space-time diagrams are plotted, representing 
the time-evolution of the vertical profiles of various hor- 
izontally averaged quantities. In addition, the location of 
unit ratio of the horizontally averaged thermal and mag- 
netic pressures, /3 = 1, is overplotted on each panel with a 
solid black line. 

Several comments can be made based on these plots. 
First, the bursty nature of the outflows and the fact that 
launching occurs approximately at the location where /3 — 1 
are both obvious from the upper left panel that shows 
the space-time diagram for the gas density. Second, the 
toroidal component of the horizontally averaged magnetic 
field {By) displays a cyclical behaviour above and below 
the midplane, a kin to the 'bu t terfly diagram' reported by 
various authors ( Grcssein 2010l : IShi et al.ll2010l : iDavis et"al] 
I2OIOI : [Simon et al. 2011b:) for stratified zero-net-flux simu- 
lations. However, in contrast to the zero-net-flux situation 
where the mean toroidal fiux arises from the disk midplane, 
it is seen here to originate at 2 ^ 2-3-ff , in agreement with 
the results reported by SI09. 

Perhaps more important are the lower left and right 
panels, where the horizontally averaged radial components 
of the magnetic field (Bx) (normalized by (i?z)) and ve- 
locity (vx) (normalized by the sound speed cq) are respec- 
tively shown (Note that, because of magnetic flux conser- 
vation, (Bz), which is used to normalize (By) in the upper 
right panel, is independent of z and t and equal to its ini- 
tial value -Bo)- They show that a significant mean radial 



magnetic field and mean radial velocity develop in the disk 
atmosphere. Their amplitudes are comparable to those of 
(Bz) and (vz), which means that both the poloidal magnetic 
field lines and the poloidal streamlines are significantly in- 
clined with respect to the vertical axis. As expected be- 
cause of the shear, the radial and toroidal components of 
the magnetic field are anti-correlated. Positive (By) corre- 
sponds to negative {B^) and vice versa. The radial veloc- 
ity also correlates positively with the radial magnetic field 
above the midplane (i.e. positive (vx) corresponds to posi- 
tive (Bx)) but correlates negatively with it below the mid- 
plane (where positive (vx) corresponds to negative {Bx)). 
This flow morphology closely r esembles that of the classi- 
cal wind solution described by iBlandford &: Pavnd (|1982l 
hereafter BP82). There are differences, though. First, in 
contrast to the steady-state BP82 solutions, the flow is 
strongly time-dependent here. Second, because of the sym- 
metries associated with the shearing-box model, outflows 
are associated with poloidal magnetic field lines inclined 
towards either positive or negative x. This is different from 
the situation usually considered in cylindrical geometry in 
which magnetic field lines are inclined away from the central 
ob ject in order t o produce a successful outflow. As noted 
bv lOgilvi^ (|2012[ ). the shearing box is unaware whether the 
central object is located in the direction of positive or nega- 
tive X, and the local physics of wind launching is symmetri- 
cal. In our shearing-box simulations, there can even be sit- 
uations where the poloidal field lines are directed towards 
positive X above the disk midplane, but towards negative x 
below the midplane (see, for example, the structure of the 
flow for 15 < t < 20 in figure [5]). Such unfamiliar flow mor- 
phologies are however permitted solutions of the equations 
in the framework of the shearing-box model ()Ogilviell2012l : 
iLesur et al.ll2012[ ) , and suggest that the outflows above and 
below the midplane are to some extent independent of each 
other. 

Because of this variety of fleld conflgurations, and since 
the flow is not in steady state, it is not straightforward to 
illustrate the analogy with the BP82 solutions more quan- 
titatively. In attempting to do so, we have isolated a single 
burst event that happens between the times t = 22 and 
t = 27 and have analysed it in more detail. Time-averaging 
the simulation data over such a short interval obviously 
means that the resulting curves will be noisy, but this is the 
price to be paid for a quantitative comparison with the the- 
oretical expectations. The structure of the time-averaged 
flow is illustrated in flgurelHl Poloidal streamlines and mag- 
netic fleld lines, spatially averaged in the y direction and 
time-averaged over the ejection event, are represented in 
the {x, z) plane. They are overplotted on the density and 
azimuthal magnetic fleld distribution respectively. Poloidal 
magnetic field lines are systematically bent towards nega- 
tive X for z > 3H and poloidal streamlines tend to align 
with those magnetic field lines, in a way that resembles 
Blandford & Payne type of solutions. In addition, a layer 
with a strong By is apparent at the base the launching re- 
gion, again consistent with this picture. In the following, 
we analysed quantitatively the properties of that flow. 

The first diagnostic we obtained is provided by the ver- 
tical mass flux. Its vertical profile is shown in figure!?] While 
small and fluctuating in the bulk of the disk, (pvz) reaches 
a constant value for \z/H\ > 3. This means the flow has 
reached a steady state in the disk atmosphere, and we can 
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Fig. 5. Space-time diagram of the logarithm of the horizontally averaged mass density (upper left panel), {By)/{Bz) 
[upper right panel), {Bx) / {Bz) [lower left panel) and {vx)/co [lower right panel) for model Diffu^H. For all panels, the 
solid line shows the location of unit ratio of horizontally averaged thermal and magnetic pressures, /3 = 1. 



write 



a(Bz 



where a is constant beyond \z\ ^ 3H. Under such con- 
ditions, it is well known that the flow is characterized 
by a number of features and invaria nts first studied in 
globa l models of disk winds (BP82; iPelletier fc Pudrit3 
Il992l) and recently ex tended to the case of the shearing 
box ( Lesur et al.ll2012l) . First, poloidal streamlines approx- 
imately align with poloidal magnetic field lines, as shown 
on figure [HI In figure [51 wc plot the vertical profiles oi 6b 
and 9y, which respectively denote the angles between the 
vertical axis and the horizontally averaged poloidal mag- 
netic and velocity fields. In those regions where the vertical 
mass flux is constant, both 9b and 9^ reach significant and 
fairly z-indcpendcnt values of order 60°. These angles are 
comparable with each other [9b ^ 9y) and larger than 30°, 
which is the critical angle in the BP82 analysis. This sup- 
ports the BP82 scenario: poloidal streamlines are approxi- 
mately aligned with poloidal magnetic fields and magneto- 
centrifugal launching is possible. Investigating whether this 
is actually the case can be examined using the Bernoulli in- 
variant, which should be constant along magnetic field lines 



(|Lesur et al.ll2012[ ): 



(5) E, 



Ber 



E 



K 



Et + Ed, 



Ep 



+ cllog[{p)) + <P - 



(B„)v, 



(6) 



where f * 



'iy) — a{By) / p and represents the shearing 
box tidal potential. The variations of the four components 
of Eser along a particular field line originating at the disk 
midplane are plotted in figure HI Their sum is also rep- 
resented using a black line, and exhibits a relatively flat 
profile beyond \z\ ^ 3H, consistent with the expectations 
of steady state theory. Acceleration of the outfiow, as seen 
by an increase of kinetic energy, occurs after the flow has 
reached a maximum in the potential energy. Its acceleration 
is essentially a result of a decrease of the potential energy 
along the field line. This demonstrates that the outflows are 
centrifugally driven, as observed for global simulations of 
laminar resistive disks. Note also that the variations of the 
thermal energy Et contribute to about one third of that 
acceleration, indicating a non-negligible effect of thermal 
pressure even above the Alfven point. This importance of 
thermal pressure is due to the plasma /3 being much larger 
in the disk atmosphere (it is of order 1) than is common in 
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Fig. 6. y-averaged poloidal streamlines {left panel) and 
poloidal magnetic field lines {right panel) in the {x,z) 
plane, time-averaged over the ejection event occurring at 
22 < t < 27, and overplotted on the density distribution 
{left panet) and azimuthal magnetic field distribution {right 
panel). 



standard studies of outflow launching (where usually /3 ^ 1 
in the launching region). 

It is well known that such disk outflows carry angular 
momentum away from the disk midplane: angular momen- 
tum is extracted from the disk by the magnetic field before 
being transferred to gas as it is accelerated (jLesur et alj 
l2012r) . This exchange is described by the following invari- 
ant that should be constant on streamlines: 



where £ = 



(7) 



- (2 — q)Q,x plays the role of specific angular 



momentum in the shearing-box approximation. Its vertical 
profile along a particular streamline is plotted in figure 1101 
Its components C and {By) /a arc also represented for com- 
parison. Beyond \z\ = 3H, f shows a flat profile, with weak 
variations that are much smaller than the variations of its 
constituents. This is again consistent with the expectations 
of BP82: disk outfl ows extract angular momentum from the 
disk. As shown by iLesur et al.l (I2012D . this leads to radial 
advection of gas and poloidal magnetic field lines. However, 
we failed to detect such an advection velocity. Averaging 
the radial mass flux over the single burst, we found 
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Fig. 7. Vertical profile of the horizontally averaged vertical 
mass flux, {pvz), obtained in model DiffuJ^H. The data are 
time-averaged over a single burst between i = 22 and t = 21 
and normalized by pqCq. 




1 1 {pvx)dtdz 
II{p)dtdz 



= 2 X 



(8) 



Fig. 8. Vertical profiles of the angles ds {red curve) and 9y 
{blue curve), measured between the vertical axis and the 
horizontally averaged poloidal magnetic and velocity fields, 
respectively. Both curves are obtained from model Diffu4H 
and are time-averaged over a single burst of ejection be- 
tween t = 22 and t = 27. 



This value is much below the turbulent velocity fluctua- 
tions, for which the Mach number is at the 10% level. This 
suggest that outflow-mediated angular momentum trans- 
port is very small, as expected given the large value of /3o 
we use, and consistent with the low altitude of the Alfven 
point. We shall return into this issue in the discussion. 

While enlightening as to the structure of the flow in 
tlie disk atmosphere, the above analysis should not hide 
the fact that tlie outflows that develop in the simulations 
are not strictly speaking BP82-type disk winds. First, one 
should not forget the non-steady nature of the flow, the 
properties of which are yet to be understood. How much of 
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Fig. 9. Vertical profiles of the quantities Ek {blue curve), 
Eb [red curve), Et {green curve), E^ [yellow curve) and 
their sum E {black) along a particular magnetic field line. 
Data have been horizontally averaged and time averaged 
between t = 22 and t = 27, i.e. over a single ejection event. 



Fig. 11. Vertical profile of the total stress tensor, normal- 
ized by the midplanc pressure, for models Diffu^H {blue 
line) and Diffu^HPmhalJ {green line). The blue and green 
dashed lines represent the corresponding values obtained in 
unstratified runs having the same parameters. 
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Fig. 10. Vertical profiles of the quantities C {blue curve), 
{By) /a {red curve) and the expected invariant f = C — 
{By) I a {green curve) along a particular magnetic field line, 
computed as in figure [S] 



that variability is an artifact of the shearing-box symme- 
tries should be investigated further. Another important dif- 
ference with classical disk wind theories is that the outflows 
presented here do not pass all the relevant critical points 
within the computational domain. It is unlikely that the 
successful complete ejection of the outflow can be demon- 
strated within the local shearing-box model because the 
depth of the gravitational potential is not represented. We 
shall return to that point in the discussion (section 15. 3p . 
To summarize, even if the comparison drawn above is ap- 
pealing, outflows driven by turbulent disks are clearly more 
complex than described by simple steady-state disk wind 
models. Working out simple models that will simultane- 



ously describe the midplane turbulent flow and the non- 
steady BP82-like winds are beyond the scope of this paper 
and will surely prove challenging in the future. 

4.2. Disk midplane: turbulent transport 

In this section, we focus on angular momentum trans- 
port. As argued above (see also section ^71\ . it is likely 
dominated by turbulent transport occurring in the bulk 
of the disk. Here, we thus present a preliminary study 
of the saturation properties of the turbulence in stratified 
shearing boxes threaded by a mean magnetic field. In ad- 
dition, we compare our results with identical simulations 
performed in the unstratified limit in order to assess the ef- 
fect of density stratification on transport properties. Given 
the large computational expense associated with such a 
survey, we defer a systematic study of the scaling with 
vertical field amplitude to a future publication and con- 
sider here only the effect of the dissipation coefficients. To 
do so, we consider two models, the parameters of which 
are respectively Re = Rm = 3000 (model Dijfu4H) and 
{Re,Rm) = (3000,1500) (model Diffu4Ha). Thus Pm = 1 
in the former and 0.5 in the latter. 

The a;y-component of the sum of tlie Reynolds and 
Maxwell stresses, normalized by the midplane pressure, 
can be taken to define a Shakura-Sunyaev a viscosity 
parameter. Its vertical profile is plotted in figure [TT] for 
both models. As reported in the literature for different 
configurations, turbulent transport increases with Pm. We 
found the vertically averaged value to be a = 2.3 x 10~^ 
when Pm ~ 1 and a = 1.1 x 10~^ when Pm = 0.5. The 
vertical profile of the stress we obtained in each case is com- 
parable to that reported in the liter ature by many authors 
using local or global simulations ([Miller &: Stond 20001: 



Hirose et all 120061: iFromang fc Nelson' '2006 



lang & ; 

20101: iDzvurkevich et al.l 1201 0: Sorat hia et all 



Flaig et all 
20101: 



Fromang et al. I 1201 ID . It is rather flat in the bulk of 
the disk, with possible local maxima at z ^ 2H before 
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Fig. 12. Same as figure 21 but for tire model having zero- 
gradient vertical boundary conditions on and By. 



dropping to low values in the corona. As shown by the 
dashed lines in figure 1111 such transport rates compare 
favourably with results obtained in unstratified boxes. We 
indeed obtained a = 4.0 x lO"'^ and 2.4 x 10"'^ when 
Pm = 1 and 0.5, respectively. These values are higher by 
only about 30% than the stratified results and display the 
same trend with magnetic Prandtl number. This means 
that unstratified sh earing-box simulations wit h net flux, 
such as reported by iLongaretti fc LesuH ()2010f ). provide a 
good first-order estimate of angular momentum transport 
in stratified disks. 



5. Discussion 

In this section, we discuss several controversial aspects 
of the present simulations, namely the influence of verti- 
cal boundary conditions on our results, the relative im- 
portance of turbulent and outflow-mediated angular mo- 
mentum transport, and possible issues associated with a 
straightforward truncation of the shearing-box vertical po- 
tential. 

5.1. Influence of the vertical boundary conditions 

Clearly, it is not satisfying to see such a strong dependence 
of the outflow mass-loss rate on numerical aspects of the 
setup. One might worry that some of this dependence is 
due to the vertical boundary conditions. In particular, our 
choice of purely vertical magnetic flelds at the top and bot- 
tom of the computational box is rather unfortunate. Indeed, 
it does not favour magneto-centrifugal acceleration as the 
latter requires magnetic field lines inclined by more than 
30° to operate. This might affect the fiow and prevent suc- 
cessful launching of disk outfiows. 

In order to test that possibility, we performed an addi- 
tional simulation with different boundary conditions on the 
magnetic field. We adopted zero-gradient boundary condi- 
tions for Bx and By at the top and bottom of the box. To 
save computational power, this simulation was performed 
without explicit dissipation coefficients and with a smaller 



resolution of 32 cells per using the same box size as for 
model DtjfuiH, namely {L^,Ly,L^) = {AH,AH,IQH). We 
obtained my, = 1.0 x 10"'^. This is very similar to the mass 
outflow rates obtained for models Diffu4H and Diffu^Ha 
that share the same box size. As an additional check, we 
plot in figure [T^ a comparison between the vertical gas ve- 
locity and the various wave velocities of the problem, as 
done for model DiffuJ^H in the top panel of figure ID The 
two plots are very similar, indicating a weak effect of the 
boundary conditions on the fiow properties (note, however 
that, as opposed to mode DijJuJ^H, the fiow velocity remains 
smaller that the sound speed everywhere). 

Thus, we have shown in this paper that our results are 
unchanged when we adopt widely used boundary condi- 
tions for the magnetic field (purely vertical field or zero gra- 
dient for its horizontal components). The good agreement 
with the results published by SI09 and SMIlO suggests that 
this insensitivity to vertical boundary conditions extends to 
more elaborate boundary conditions. This is strange as we 
have shown that the results depend on vertical box size 
(section I3.4.2P which indicates that the flow knows about 
the existence of an upper boundary. Its seems that the ar- 
tifacts uncovered in the present paper are due to the very 
existence of a boundary rather that the nature of the con- 
ditions imposed there. Alternatively, this could indicate a 
breakdown of the shearing-box model itself. Future investi- 
gations along with comparison with global simulations are 
needed to clarify this issue. 

5.2. Vertical vs. turbulent transport of angular momentum 

The paradigm that emerges out of the analysis performed in 
this paper is that of a turbulent disk surrounded by tran- 
sient and presumably local magneto-centrifugally acceler- 
ated BP82-like outflows. As demonstrated in flgurefTUl such 
outflows carry angular momentum away from the bulk of 
the disk, the amount of which might be important for its 
evolution. In this section, we speculate on the relative im- 
portance of the angular momentum transported vertically 
by the outflow and radially by the turbulence. 

Given the symmetries of the shearing box, there is some 
ambiguity in determining that ratio. One possibility we ex- 
plored is to separate the flow into mean fields and fluctu- 
ations and to assign the transport mediated by the former 
to the outflow and that due to the latter to the turbulence. 
However, we found significant radial transport due to the 
mean magnetic field, the role of which remai ns unclear. We 
therefore followed the approach sketched by I Wardl^ (|2007ft 
for a radially structured accretion disk. In cylindrical co- 
ordinates {R,(j),Z), we first write the horizontal velocities 
Vfi and as the sum of an azimuthally averaged part and 
remaining fluctuations: 

VB. = {vb) + 5vR , (9) 
V4, = {v<i>) + Sv^ = Rfl + 6v^ . (10) 

Angular momentum conservation in a steady-state disk is 
then written 



R{pvR)^rjr - 7m (^'^^-a) + {R'Tz,) , (11) 
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with the turbulent stress ten sors Tr^ and T^^ expressed as 
(|Balbus fc Papaloizoul[T999[l : 



Tr4 = {BrB^ - pSvrSv^) 



(12) 
(13) 



Wc now assume that there exist disk surfaces located at 
Z = ±Zs{R) above and below which mass and angular 
momentum are outflowing. We define the mass accretion 
rate in the disk as 



M 



{pvR)dZ (14) 
and integrate Eq. ([TT|) between —Z^ and +Zs. We find 



2-kR dR 



ruR + mz 



(15) 



where ttir and mz respectively denote the disk and the 
outfiow contribution to the angular momentum transport 
and are given by 



mz = -R[Tz4,]tt 



Tr^cIZ 



(16) 
(17) 



The relative importance of turbulence and outflowing mate- 
rial to angular momentum transport is given by the ratio of 
TTiR and mz- In the shearing box, however, m^j artiflcially 
vanishes because of symmetries of the model. It is possible, 
though, to use the shearing box results to make such a eval- 
uation for real disks, provided reasonable assumptions are 
made about the disk structure. Let us first assume that the 
vertical profile of the stress, when normalized by the mid- 
plane pressure Pmid, is a "universal" function oi Z = Z/H 
only, and is independent of R (which amounts to saying 
that a is radially constant). Then all the radial dependence 
in niR can be made explicit provided we give ourselves ra- 
dial scalings for the surface density S cx i?P and the sound 



speed 



R'^ (we assume a locally isothermal equation 



of state, i.e. that cq a function of R only, which means the 
vertical density profile is Gaussian). After some algebra, we 
obtain 

TOi? - (p + <z + 2)a^-=^, (18) 



27r 

where we have defined the radius-independent variable Wr 
as 

/ aR{Z)dZ ^ - f / ^dZ, (19) 



P, 



■mid 



and the second equality serves as a definition of aR{Z). 
The interest of the last pair of equations is that cxr can be 
derived from shearing-box simulations such as presented in 
the present paper. Similarly, rhz can be written as 



mz 



V2^' 



(20) 



where az = —Tz^/Pmid, so that the ratio mR/mz is given 
by _ 

b+^+2)rf)^^. (21) 
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Fig. 13. Vertical profile of ur [blue curve) and az {red 
curve) , time averaged using the data of model Diffu4-H over 
the interval 22 < i < 27. Vertical arrows marks the location 
of unit ratio between horizontally averaged thermal and 
magnetic pressures. 



Note the absolute value in the expression above which ac- 
counts for the facts that outflow-mediated angular momen- 
tum transport can be directed either way in the shearing 
box. We are interested here only in its amplitude. 

Using model Diffu4H, we measured ur^ and az- We 
focus on the same burst that we studied in section 14.11 
and thus averaged the simulation data between t = 22 
and t ~ 27. The vertical profiles we obtained by doing so 
are plotted in figure [13] along with vertical arrows showing 
the ejection region calculated as in fig. [31 At first glance, 
fig. [T51 shows than aR is larger than az- This is confirmed 
by precise measurements, which give or = 1.0 x 10~^ and 
[az]1lzl = -4.4 X 10"^. Note that we have taken Zs = 2-5H 
in making that estimate. Taking p + 17 + 2 to be of order 
unity in Eq. (|2T|) . this gives 



mR 



24 

\mz\ \R 



(22) 



This means that turbulent transport dominates outflow- 
mediated transport as long as H/R > 0.04. We stress, 
however, that outflows are more difficult to launch in thin- 
ner disks because the escape Mach number is larger in 
that case. Indeed, v^sc ~ V^v^, which translates into 
Vesc/co = ^/2{R/H)- Thus, in thin disks where we can 
imagine outflow-mediated transport to dominate, its ver- 
tical Mach number should exceed ^ 35 for successful ejec- 
tion. As demonstrated in both panels of figure [H this is far 
from being the case in the simulations we have been able 
to run, in which the vertical Mach number reaches at most 
values of a few. Thus the question of whether such outflows 
will eventually escape the potential well or fall back on the 
disk is still an open issue. 

5.3. Smoothing the shearing-box potential? 

It is obvious that the shearing-box model suffers from a 
number of artifacts that might affect the results presented 
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here. Maybe the most important of them is the singular 
nature of the tidal potential used in our simulations. In 
standard shearing-box simulations, the potential is given 
as 

4,^, = -^n^x^ + ^n^z\ (23) 

The second term shows that the potential becomes deeper 
and deeper as the vertical box size is increased. This means 
it requires larger and larger energy for gas to climb out of 
that potential, while in reality the gravitational potential 
has a finite depth 



(24) 



The second term in Eq. (j23p is simply the second-order 
Taylor expansion of the last equation with respect to z. 

SMIlO investigated this issue by directly employing the 
potential given by Eq. ([M)) . They obtain results that are 
difficult to interpret as liiw display variations by factors of 
2-3 and show no systematic trend with vertical box size. 
Here, we caution that a more rigorous and self-consistent 
approach is to include the next order in the expansion of the 
tidal potential. Up to the first terms that include a modifi- 
cation of its vertical profile, such an expansion is written 



-1, ^ 2 

X XZ 

2 



(25) 



where we have noted x = x/H and z ~ z/H and use the 
relation H = co/fl. 

The expression above shows that the first additional 
terms in the expansion involves x as well as z. This means 
that smoothing the potential should be done by simulta- 
neously altering the radial component of the tidal force. 
We note that in such a situation, curvature terms should 
also be accounted for, and the shearing-periodic boundary 
conditions would be incompatible with the governing equa- 
tions. The first nonzero term involving z that appears in the 
potential expansion scales like xz'^. In practice, this term 
change the equilibrium shear itself, making it z-dependent. 
For a barotropic gas as considered here, this violates the re - 
sults that rotation should be on cylinders (lTassoullll978l) . 
unless the equilibrium profile of the density depends on x as 
well as on z. This is a significant modification to the stan- 
dard shearing box model. In practice, smoothing can only 
be achieved with higher order terms that are in (H/R)^, 
among which a term in — z* that results in a finite poten- 
tial barrier out of which to climb, albeit with an incorrect 
value. These complications were not addressed in the SMIlO 
analysis, which makes the interpretation of their results un- 
clear. 

Clearly, the discussion above shows that changing the 
vertical profile of the tidal potential is a risky avenue, much 
beyond the scope of this paper. It has to be done in a way 
that is consistent with all the approximations involved in 
the shearing box. The extension of the shearing-box model 
in that direction is promising for the analysis of outflow- 
related issues and should be the subject of future work. 



6. Conclusions and perspectives 

In this paper, we have analysed in detail the structure of 
the flow in an accretion disk threaded by a weak vertical 



magnetic field. In agreement with previous results (SI09, 
SMIlO), we find that a strong outflow develops, with mass- 
loss rates that deplete the gas content of the computa- 
tional box within a hundred dynamical times. However, we 
found these mass-loss rates to depend strongly on the nu- 
merical setup, thus preventing a reliable estimate of such 
outflow rates from being made in shearing-box simula- 
tions. Nevertheless, we found that the flow properties in 
the disk atmosphere exhibit robust features: outflows are 
magneto-centrifugally accelerated, yet thermal pressure has 
a non negligible effect. We recover the same invariants that 
charact erized classical theories o f steady-state disk winds 
fBP82: lPelletier fc Pudr"it3ll992[ ). We stress however, that 
there are a number of important differences with such the- 
ories: first, the outflows produced in the simulations are 
strongly time-dependent, which means steady state models 
should be used with care. Second, the vertical outflow veloc- 
ity is never seen to exceed the fast magnetosonic speed. In 
addition, the positions of the slow magnetosonic point and 
the Alfven point are seen to depend on box size. This high- 
lights the fact that information can come back down the 
outflow, although we have not been able to flnd any signifi- 
cant modifications of its properties by changing the nature 
of the vertical boundary conditions. Finally, and maybe 
most importantly, we found that the dynamical effect of 
these outflows on accretion disks remains to be demon- 
strated: we investigated the ratio of the angular momen- 
tum extraction by such outflows and its radial transport 
by turbulence. These quantities might become comparable 
for thin disks for which the escape velocity is larger by more 
than one order of magnitude than the velocities observed in 
the simulations. Thus, the ultimate fate of these outflows, 
either successful launching or falling back onto the disk, 
remains unclear. 

These results should be contrasted with a case of 
stronge r , MR I stable magnetic fleld. In that situation, 
lOgilvid (|2012[ ) has argued, with reference to previous work, 
that the local approximation can be used to determine the 
mass-loss rate per unit area in magnetized outflows that 
are accelerated along incli ned poloidal magnetic field lines 
(jBlandford fc Pavnd 11983) . In turbulent disks as studied 
in the present paper, this result does not hold. The oc- 
currence of transient, and presumably local, outflows ap- 
pears to be a natural consequence of the nonlinear devel- 
opment of the MRI in a stratified disk in the presence of 
a net vertical magnetic field. This behaviour is different 
from that in incompressible disks, where no vertical motion 
arises from nonlinear channel modes bec ause the vertical 
force s are cancelled by pressure gradients (jGoodman fc Xul 
I1994D . Nevertheless, the connection between these transient 
and local outflows and the large-scale jets observed from 
astrophysical discs is still unclear. Traditional models for 
large-scale outflows have strong, ordered magnetic fields 
and Alfven points high above the surface of the disk, rather 
different from the situation in our simulations. 

In the light of our results, it appears that the shearing- 
box model is not well adapted to answer these questions, 
unless new developments arc undertaken that go much be- 
yond its original formulation. Both quantitative mass-loss 
rates and the ultimate fate of these outflows may have to 
be determined using global numerical simulations. This is 
a formidable task, the computational cost of which is at 
the limit of currently available facilities. But perhaps the 
good news of our work is that moderate resolution and run- 
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time are sufficient to capture the properties of disk outflows. 
This gives hope that the questions of the fate, properties 
and consequences of outflows for accretion disk dynamics 
can be answered in the next few years. 
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Appendix A: Linear modes in viscous and resistive 
stratified disks 

Consider an MRI channel flow of the following form: 



u 



ch 



B 



ch 



Hne''*[iu^{z) + kuy{z)], (A.l) 
Boe''[iB,{z) + kBy{z)], (A.2) 



where H is the disk semithickness, Bq is the amplitude of 
the background vertical field, s is the growth rate, and u^, 
Uy, Bx and By are dimensionless functions to be deter- 
mined. Note that in the limit of high midplane /3 (anelastic 
approximation) this flow is both a linear and nonlinear so- 
lution to the equations of non-ideal MHD in the vertically 
stratified shearing box (cf. iLatter et al.ll2010l) . 

We denote the background density as p = po h{z) where 
/i is a dimensionless function, and then scale time by 1/VL 
and space by H . The four equations to be solved are then 



2ii,, 



S Uy 

sB:c 

sB„ 



2 1, 11 

P h Reh 

Rm 



■ B^ 



Rm ^ ' 



(A.3) 
(A.4) 
(A.5) 
(A.6) 



where a prime indicates differentiation with respect to z. 
This set of equations can be solved via a pseudo-spectral 
method using Whittaker cardinal functions, as outlined in 
Latter et al. (2010). 

Examples of these calculations are plotted in Fig. lA.ll 
Note in these profiles how the combined action of magnetic 
diffusion and viscosity dramatically alters the mode mor- 
phology even for relatively large Re and Rm. In the ideal 
case, the modes are concentrated near the miplane, where 
they also exhibit their shortest scales. These short scales, 
however, are vulnerable to even a small amount of dissipa- 
tion, and when present the modes as a consequence tend to 
localise near z — zL2-3H where their characteristic scale is 
longer. 
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Fig. A.l. Vertical profile of the radial magnetic field eigen- 
mode for three different sets of dissipation coefficients. 
From left to right, Rm = Re = 10^, Rm = Re = 10^ 
and Rm = Re = 3000. 



